Mixed Model for Photosynthesis in Ulva lactuca
Angela Richards DonĂ
5/17/2022
Run 5.6 PHOTOSYNTHESIS Analysis, Script Chunks, and Plots
This is the analysis of the final run of the Ulva lactuca salinity and nutrient experiments conducted on the lanai in St. John 616. These experiments incorporated four paired salinity and nutrient levels with three temperature levels. Each run produced an n = 2 and was repeated initially 8 times for a total of n = 16. Data gaps were identified and filled by early March 2022. This output reflects all data totaling five treatments for Ulva lactuca.
Packages loaded:
library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(ggplot2)
library(ggpubr)Load and prepare the dataset
Open the output dataset generated by the ps_script_clean_to_ek_alpha.R script in the phytotools_alpha_ek project
run5_6_photosyn_data <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/run5-6_ek_alpha.csv")Assign run as a factor
run5_6_photosyn_data$Run <- as.factor(run5_6_photosyn_data$Run)Assign temperature as a factor
run5_6_photosyn_data$Temperature <- as.factor(run5_6_photosyn_data$Temp...C.)Assign treatment as characters from integers then to factors
run5_6_photosyn_data$Treatment <- as.factor(as.character(run5_6_photosyn_data$Treatment))Subset the data and toggle between the species for output. Use Day 9 for final analysis ONLY
ulva <- subset(run5_6_photosyn_data, Species == "ul" & RLC.Day == 9)Add a column for growth rate from growth rate dataset to the alredy subsetted ulva data frame
growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/run5-6_growth_all_042922.csv")gr_ulva <- subset(growth_rate, Species == "Ul")
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)rETRmax – Run the model
Run model for rETRmax with two fixed effect variables and three random effects variables
run5_6_photosyn_model_noint <- lmer(formula = rETRmax ~ Treatment + Temperature
+ (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order),
data = ulva)rETRmax – Make a histogram and residual plots of the data
hist(ulva$rETRmax, main = paste("Ulva lactuca rETRmax"), col = "darkolivegreen2", labels = TRUE)
plot(resid(run5_6_photosyn_model_noint) ~ fitted(run5_6_photosyn_model_noint))
qqnorm(resid(run5_6_photosyn_model_noint))
qqline(resid(run5_6_photosyn_model_noint))
rETRmax – Check the performance of the model
performance::check_model(run5_6_photosyn_model_noint)
These outputs show the model is acceptable
rETRmax – Check r2 for model fit and print the model statistics summary
r.squaredGLMM(run5_6_photosyn_model_noint)## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.## R2m R2c
## [1,] 0.5882889 0.6835032summary(run5_6_photosyn_model_noint)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rETRmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
## (1 | RLC.Order)
## Data: ulva
##
## REML criterion at convergence: 1827.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8253 -0.6047 0.0845 0.5200 3.4082
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 22.155 4.707
## RLC.Order (Intercept) 5.044 2.246
## Run (Intercept) 6.305 2.511
## Residual 111.367 10.553
## Number of obs: 240, groups: Plant.ID, 95; RLC.Order, 6; Run, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 34.819 3.437 3.609 10.130 0.000887 ***
## Treatment1 20.641 3.669 3.083 5.626 0.010318 *
## Treatment2 21.727 3.669 3.083 5.922 0.008906 **
## Treatment3 37.848 3.669 3.083 10.317 0.001724 **
## Treatment4 39.172 3.669 3.083 10.677 0.001555 **
## Temperature27 -4.430 2.504 30.316 -1.770 0.086835 .
## Temperature30 -2.157 2.293 65.525 -0.941 0.350296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.721
## Treatment2 -0.721 0.828
## Treatment3 -0.721 0.828 0.828
## Treatment4 -0.721 0.828 0.828 0.828
## Temperatr27 -0.348 0.003 0.003 0.003 0.003
## Temperatr30 -0.338 0.000 0.000 0.000 0.000 0.475rETRmax – Run an ANOVA and pairwise comparisons based on ANOVA results Results for temperature are non-significant, thus no pairwise comparisons for this
anova(run5_6_photosyn_model_noint, type = c("III"), ddf = "Satterthwaite")ulva_photosyn_model_aov <- aov(rETRmax ~ Treatment + Temperature, data = ulva)
TukeyHSD(ulva_photosyn_model_aov, "Treatment", ordered = FALSE)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = rETRmax ~ Treatment + Temperature, data = ulva)
##
## $Treatment
## diff lwr upr p adj
## 1-0 20.610833 13.932211 27.289455 0.0000000
## 2-0 21.696250 15.017628 28.374872 0.0000000
## 3-0 37.818125 31.139503 44.496747 0.0000000
## 4-0 39.141667 32.463045 45.820289 0.0000000
## 2-1 1.085417 -5.593205 7.764039 0.9917098
## 3-1 17.207292 10.528670 23.885914 0.0000000
## 4-1 18.530833 11.852211 25.209455 0.0000000
## 3-2 16.121875 9.443253 22.800497 0.0000000
## 4-2 17.445417 10.766795 24.124039 0.0000000
## 4-3 1.323542 -5.355080 8.002164 0.9824820rETRmax – Plot the results
plot(allEffects(run5_6_photosyn_model_noint))
Plot a regression between the photosynthetic independent variables of interest and growth rate
ulva_growth_etr_graph <- ggplot(ulva, aes(x=rETRmax, y=growth_rate)) + geom_point() +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Ulva lactuca rETRmax vs Growth Rate", x = "rETRmax (μmols electrons m-2 s-1)",
y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_graph## `geom_smooth()` using formula 'y ~ x'
Temperature did not have a significant effect on the outcome for rETRmax but it did for Ek and alpha. alpha histogram is NOT Normal, probably not a good fit for this analysis
Ek – Run the model
Run model for minimum saturating irradiance (Ek) with two fixed effect variables and three random effects variables
run5_6_photosyn_model_noint <- lmer(formula = ek.1 ~ Treatment + Temperature
+ (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order),
data = ulva)Ek – Make a histogram and residual plots of the data for ulva
hist(ulva$ek.1, main = paste("Ulva lactuca Ek"), col = "darkolivegreen3", labels = TRUE)
plot(resid(run5_6_photosyn_model_noint) ~ fitted(run5_6_photosyn_model_noint))
qqnorm(resid(run5_6_photosyn_model_noint))
qqline(resid(run5_6_photosyn_model_noint))
Ek – Check the performance of the model
performance::check_model(run5_6_photosyn_model_noint)
Ek – Check r2 for model fit and print the model statistics summary
r.squaredGLMM(run5_6_photosyn_model_noint)## R2m R2c
## [1,] 0.5867698 0.7012116summary(run5_6_photosyn_model_noint)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ek.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
## (1 | RLC.Order)
## Data: ulva
##
## REML criterion at convergence: 2102.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.04017 -0.53517 -0.01583 0.51382 2.98662
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 66.38 8.147
## RLC.Order (Intercept) 10.92 3.304
## Run (Intercept) 61.94 7.870
## Residual 363.52 19.066
## Number of obs: 240, groups: Plant.ID, 95; RLC.Order, 6; Run, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 16.796 8.856 3.225 1.897 0.14770
## Treatment1 30.241 9.766 3.068 3.096 0.05187 .
## Treatment2 40.947 9.766 3.068 4.193 0.02368 *
## Treatment3 68.624 9.766 3.068 7.027 0.00550 **
## Treatment4 71.137 9.766 3.068 7.284 0.00495 **
## Temperature27 -9.826 4.311 29.565 -2.279 0.03004 *
## Temperature30 -8.349 4.020 68.996 -2.077 0.04154 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.820
## Treatment2 -0.820 0.921
## Treatment3 -0.820 0.921 0.921
## Treatment4 -0.820 0.921 0.921 0.921
## Temperatr27 -0.235 0.002 0.002 0.002 0.002
## Temperatr30 -0.229 0.000 0.000 0.000 0.000 0.480Ek – Run an ANOVA and pairwise comparisons based on ANOVA results
anova(run5_6_photosyn_model_noint, type = c("III"), ddf = "Satterthwaite")ulva_photosyn_model_aov <- aov(ek.1 ~ Treatment + Temperature, data = ulva)
TukeyHSD(ulva_photosyn_model_aov, "Treatment", ordered = FALSE)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ek.1 ~ Treatment + Temperature, data = ulva)
##
## $Treatment
## diff lwr upr p adj
## 1-0 30.19583 17.895251 42.49642 0.0000000
## 2-0 40.90208 28.601501 53.20267 0.0000000
## 3-0 68.57917 56.278584 80.87975 0.0000000
## 4-0 71.09167 58.791084 83.39225 0.0000000
## 2-1 10.70625 -1.594333 23.00683 0.1208714
## 3-1 38.38333 26.082751 50.68392 0.0000000
## 4-1 40.89583 28.595251 53.19642 0.0000000
## 3-2 27.67708 15.376501 39.97767 0.0000000
## 4-2 30.18958 17.889001 42.49017 0.0000000
## 4-3 2.51250 -9.788083 14.81308 0.9803988TukeyHSD(ulva_photosyn_model_aov, "Temperature", ordered = FALSE)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ek.1 ~ Treatment + Temperature, data = ulva)
##
## $Temperature
## diff lwr upr p adj
## 27-20 -9.9375 -18.11226 -1.7627403 0.0125118
## 30-20 -9.0075 -17.18226 -0.8327403 0.0267837
## 30-27 0.9300 -7.24476 9.1047597 0.9610888Ek – Plot the results
plot(allEffects(run5_6_photosyn_model_noint))
ulva_growth_ek_graph <- ggplot(ulva, aes(x=ek.1, y=growth_rate)) + geom_point() +
geom_smooth(method = "lm", col = "black") + theme_bw() +
labs(title = "Ulva lactuca Ek vs Growth Rate", x = "Ek (μmols photons m-2 s-1)", y = "growth rate (%)") +
stat_regline_equation() + stat_cor(label.y = 160)
ulva_growth_ek_graph## `geom_smooth()` using formula 'y ~ x'
alpha – Run the model
Run model for alpha withtwo fixed effect variables and three random effects variables
run5_6_photosyn_model_noint <- lmer(formula = alpha.1 ~ Treatment + Temperature
+ (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order),
data = ulva)## boundary (singular) fit: see help('isSingular')alpha – Make a histogram and residual plots of the data for ulva
hist(ulva$alpha.1, main = paste("Ulva lactuca alpha"), col = "darkolivegreen4", labels = TRUE)
plot(resid(run5_6_photosyn_model_noint) ~ fitted(run5_6_photosyn_model_noint))
qqnorm(resid(run5_6_photosyn_model_noint))
qqline(resid(run5_6_photosyn_model_noint))
alpha – Check the performance of the model
performance::check_model(run5_6_photosyn_model_noint)
alpha – Check r2 for model fit and print the model statistics summary
r.squaredGLMM(run5_6_photosyn_model_noint)## R2m R2c
## [1,] 0.4524773 0.5264424summary(run5_6_photosyn_model_noint)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: alpha.1 ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +
## (1 | RLC.Order)
## Data: ulva
##
## REML criterion at convergence: 453
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5531 -0.5754 -0.1204 0.3199 3.6187
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plant.ID (Intercept) 0.01181 0.1087
## RLC.Order (Intercept) 0.00000 0.0000
## Run (Intercept) 0.04206 0.2051
## Residual 0.34486 0.5872
## Number of obs: 240, groups: Plant.ID, 95; RLC.Order, 6; Run, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.15229 0.22956 3.28704 9.376 0.00175 **
## Treatment1 -0.84972 0.25969 3.44592 -3.272 0.03822 *
## Treatment2 -1.11836 0.25969 3.44592 -4.306 0.01731 *
## Treatment3 -1.48303 0.25969 3.44592 -5.711 0.00723 **
## Treatment4 -1.52872 0.25969 3.44592 -5.887 0.00656 **
## Temperature27 0.31167 0.09830 42.27696 3.171 0.00283 **
## Temperature30 0.25233 0.09804 46.29620 2.574 0.01332 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmn1 Trtmn2 Trtmn3 Trtmn4 Tmpr27
## Treatment1 -0.830
## Treatment2 -0.830 0.893
## Treatment3 -0.830 0.893 0.893
## Treatment4 -0.830 0.893 0.893 0.893
## Temperatr27 -0.214 0.001 0.001 0.001 0.001
## Temperatr30 -0.214 0.000 0.000 0.000 0.000 0.499
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')alpha – Run an ANOVA and pairwise comparisons based on ANOVA results
anova(run5_6_photosyn_model_noint, type = c("III"), ddf = "Satterthwaite")ulva_photosyn_model_aov <- aov(alpha.1 ~ Treatment + Temperature, data = ulva)
TukeyHSD(ulva_photosyn_model_aov, "Treatment", ordered = FALSE)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = alpha.1 ~ Treatment + Temperature, data = ulva)
##
## $Treatment
## diff lwr upr p adj
## 1-0 -0.8506250 -1.1973752 -0.50387476 0.0000000
## 2-0 -1.1192708 -1.4660211 -0.77252059 0.0000000
## 3-0 -1.4839375 -1.8306877 -1.13718726 0.0000000
## 4-0 -1.5296250 -1.8763752 -1.18287476 0.0000000
## 2-1 -0.2686458 -0.6153961 0.07810441 0.2109076
## 3-1 -0.6333125 -0.9800627 -0.28656226 0.0000101
## 4-1 -0.6790000 -1.0257502 -0.33224976 0.0000018
## 3-2 -0.3646667 -0.7114169 -0.01791643 0.0338202
## 4-2 -0.4103542 -0.7571044 -0.06360393 0.0113476
## 4-3 -0.0456875 -0.3924377 0.30106274 0.9963020TukeyHSD(ulva_photosyn_model_aov, "Temperature", ordered = FALSE)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = alpha.1 ~ Treatment + Temperature, data = ulva)
##
## $Temperature
## diff lwr upr p adj
## 27-20 0.303925 0.07348064 0.5343694 0.0059250
## 30-20 0.252950 0.02250564 0.4833944 0.0275126
## 30-27 -0.050975 -0.28141936 0.1794694 0.8607823alpha – Plot the results
plot(allEffects(run5_6_photosyn_model_noint))
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.